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Abstract 

Spray jet in crossflow has been a subject of research because of its wide application in 
systems involving pollutant dispersion, jet mixing in the dilution zone of combustors, and 
fuel injection strategies. The focus of this work is to investigate dispersion of a 
2-dimensional atomized spray jet into a 2-dimensional crossflow. A quick computational 
method is developed using available software. The spreadsheet can be used for any 2D 
droplet trajectory problem where the drop is injected into the free stream eventually 
coming to the free stream conditions. 

During the transverse injection of a spray into high velocity airflow, the droplets (carried 
along and deflected by a gaseous stream of co-flowing air) are subjected to forces that 
affect their motion in the flow field. Based on the Newton’s Second Law of motion, four 
ordinary differential equations were used. These equations were then solved by a 
4 th -order Runge-Kutta method using Excel software. 

Visual basic programming and Excel macrocode to produce the data facilitate Excel 
software to plot graphs describing the droplet’s motion in the flow field. This program 
computes and plots the data sequentially without forcing users to open other types of 
plotting programs. A user’s manual on how to use the program is also included in this 
report. 


Nomenclature 

2 

Ad Projected area of the droplet, K r 

Cd Drag coefficient, see equation (8) 

F Force 
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g Gravitational acceleration constant, 9.8 m/ 5 2 
n n th iteration 

Re Reynolds number, see equation (9) 
r Spherical radius of the droplet 
u g Velocity of the cross stream (air) in the x-direction 

Ud Velocity of the droplet in the x-direction 

— > 

U Velocity 

— > 

XJ R Relative velocity between the droplet and the gas stream 

Ur = {u d - u g)i + ( w d~ w g )k 
4 

Vd Droplet Volume, — ;rr 3 

w g Velocity of the cross stream (air) in the z-direction 

Wd Velocity of the droplet in the z-direction 

x horizontal direction ± equals to the right 
z vertical direction ± equals up 

p g Density of the cross stream (air) 

Pd Density of the droplet 

jig Viscosity constant of the gaseous fluid 

At time step-size 

SUBSCRIPTS: 
d Droplet 

g Gaseous phase 


Introduction and Governing Equations 

A liquid spray injected into a gaseous crossflow is important because of its wide 
application in systems involving two phase mixing and in combustion requiring quick 
mixing and reduction of pollutants, for jet mixing in the dilution zone of combustors, and 
for determining fuel injection strategies. It is important to be able to compute this 
flow to optimize the mixing strategy. 

This work is mainly focused on producing a quick computational method for determining 
spray penetration. With this spreadsheet, one can investigate the dispersion of an air-blast 
atomized spray jet into a crossflow, see Figure 1. During the transverse injection of a 
spray into high velocity airflow, the droplets (carried along in the gaseous stream of co- 
flowing air) are subjected to forces that affect their motion in the flow field (see 
Figure 2). 
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Figure 1 (from NASA/CR — 2000-210467, reference 1) 



The trajectories of the droplets can be tracked by applying a Lagrangian-based analysis to 
the droplets. Since no evaporation is assumed, the code can be used for solid particles as 
well. The momentum equations for a droplet can be obtained by equating the droplet 
motion to: (1) The viscosity and pressure-related drag forces, (2) The pressure gradient 
and viscous forces related to the fluid surrounding the droplet, (3) The inertia of the 
virtual mass, induced when the particle acceleration affects the fluid mass acceleration, 
and (4) The Basset force, which takes into account the acceleration history of the droplet. 


Based on these principles along with the following assumptions: 


(1) The droplets are spherical, 

(2) No droplets breakup occurs, 
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(3) Vaporization is not considered and is assumed negligible, and 

(4) Lift, virtual mass, and Basset forces are neglected. 

(5) Chemical reaction is not included, 

droplet trajectory and velocity with respect to time can be calculated. These assumptions 
reduce the droplet momentum equation to include only the effects of the drag and body 
forces. The general momentum equations for a single droplet injected along the positive 
x-direction, transversely into a downward-flowing air stream in the positive z-direction, 
as shown in Figure 2, is described by 

Fd ~~ ^drag "*■ ^body 0 ) 


where the net force F d that drives the droplet motion is balanced by the drag force 

opposing its motion, and the field forces acting on the droplet. The aerodynamic drag 
force is given by 


^ drag 2 P R 


u 


R 


4 iC D 


( 2 ) 


where p„ is the air density, and Ad and C D , the projected area and the drag coefficient of 
the droplet, respectively. The relative velocity between the droplet and the crossflow has 
a magnitude of Ur (see Figure 2). The subscript “d” refers to the droplet and “g” the 
crossflow air. 

The body force, resulting from an equivalent volume of air that buoys the droplet, 
includes the gravitational and buoyancy forces. It is given by 

F b ody=(p g -pd) V d g ( 3 ) 

which says that the body force is equal to the product of relative droplet and air density 
( pd - p g ), the droplet volume V ( /, and the gravitational acceleration g. 

Substituting equation 2 and equation 3 to equation 1 yields: 



Pd^d 


dW d 

dt 



-w) 


U t 


Ad Cd 


+ 


fa 




( 5 ) 

( 6 ) 
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(7) 


dz 

dt 


= w„ 


The drag coefficient of the droplet depends on the droplet Reynolds number and is given 

by 




24 


Re 


d L 


1 , 1 r, 2/3 

! + - Re </ 

6 


0.424 


R& d < 1000 


Re d >1000 


( 8 ) 


where Re,/ is the droplet Reynolds number and is defined as follows 


Re^ = 


2P, 


U 


R 


A, 


(9) 


in which r d is the droplet radius and jig is the gas (air) viscosity. 


Numerical Method 

Four ordinary differential equations are to be solved, namely, equations 4, 5, 6 and 7, for 
the four dependent variables uj, Wd, x and z. The droplet trajectory is defined by the set of 
x and z values. A 4 th -Order Runge-Kutta explicit method 1 was used to solve these 
equations. The Runge-Kutta explicit method is an ideal numerical scheme for solving 
ordinary differential equations using Excel software. It is a self-starting method with 
good stability characteristic. The step-size can be changed as desired without any 
complications for higher-order schemes. For a set of two coupled equations, such as, 


dx 

— = f(x,z,t) 
dt 

(10a) 

dz 

— = g(x,z,t ) 
dt 

(10b) 


the 4 th -order Runge-Kutta method reads (subscript n stands for the n th time step); k and l 
are unknown constants. 


1 


Nakamura, Shoichiro, Applied Numerical Method With Software , Englewood Cliffs, New Jersey: Prentice Hal, 1991. 
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(11a) 


X n+\ ~ X n + t (^ 1 +2 ^2 + 2 ^3 + ^4 ) 

6 

z »+i = z » + “ (h +2 h + + h ) 

6 


where 


and 


*i =A t- f(x n ,z n ,t n ) 

, a rs K U A*. 
k 2 = Af • /(*„ + y , z „ + yf„ + y) 

7 A /., K L At. 

k 3 =At- f(x n + -,z n +-,t n +—) 
k 4 =At- f(x n +k 3 ,z n +l 3 ,t n +At) 

l l =At-g(x n ,z n ,t n ) 

, . , k, l At. 

h ~ At ’ &(*« + y > z „ + — A + — ) 

7 * / K L At. 

h ~ At ' &( x « + y > z « + — A + y) 

h = At-g{x n +k 3 ,z n +l 3 ,t n +At) 


(lib) 


(12a) 


(12b) 


Only every nth cycle (as specified by the user) is saved for plotting. This greatly saves on 
storage and increases the speed of post processing. We have chosen to enter the data in SI 
units in the unlocked cells. The required conversions are done in the locked cells. When 
the user becomes familiar with the spreadsheet, the spreadsheet can be unlocked, with 
password NASA, and the user can adapt the spreadsheet as required. 

The predictions were compared with the data of reference 1 with good results, see 
Figures 3 and 4. The experimental data from Figure 9.4 in ref. 1 for cross flow jets 
without an air blast assist was used. Jet to crossflow momentum flux ratio is used in this 
study to determine the depth of the droplet penetration. Momentum flux ratio for single 
phase jet is given by 


pU : 


<h = 


jet 


pU : 


crossflow 
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1 atm 


3 atm 


5 atm 









Figure 3 — Experimental Comparison 
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Or Fortran 

I ' - ► Excel 


-0.5 


Droplet 




radius 

= 50 pm * * 

-1 


viscosity 

= 0.000018 Kg/m-s 



g 

= 9.8 m/s 2 



pd 

= 822 Kg/m 3 

-1.5 


P< 

= 1.205 Kg/m 3 

i 


At 

= 0.000001 sec 

-2 


Ug 

= 0m/s 



Wg 

= -132 m/s 

-2.5 


Initial 




Droplet Position (x, z) = (0, 0) m 



Droplet Velocity(Ud,Wd) = (0.3,0) m/s ( 

-3 



‘ 

-3.5 


i . i i i 



< 

3 

0.0002 

0.0004 0.0006 0.0008 


V2 


-0.5 7 

-1 j 

-1.5 7 

•D 

* -2 7 
-2.5 7 
-3 7 
-3.5 1 


Droplet 

radius = 50 pm 

viscosity =0.000018 Kg/m-s 
g =9.8 m/s 2 

pd =400 Kg/m 3 

p« =1.205 Kg/m 3 

At =0.000001 sec 

Ug =0m/s 

Wg =-132 m/s 

Initial 

Droplet Position (x, z) = (0, 0) m 
Droplet Velocity(Ud, Wd) = (0.3, 0) m/s 


0.0002 0.0003 

V2 


Fortran 

Excel 



i i l i 

0.0004 


1 

l 

l 

I 

I 

I 

I 

-U 


-0.5 

-1 

-1.5 

CD 

> -2 
-2.5 
-3 
-3.5 


- Droplet 


radius = 25 pm 
viscosity =0.00001 8 Kg/m-s 


g 

pa 

P* 

At 

Ug 

Wg 


= 9.8 m/s 2 
= 400 Kg/m 3 
= 1.205 Kg/m 3 
= 0.000001 sec 
= 0m/s 
= -132 m/s 


Fortran 

Excel 


Initial 

Droplet Position (x, z) = (0, 0) m 
Droplet Velocity(Ud, Wd) = (0.3, 0) m/s 


V2 


Figure 4 — Droplet Trajectory Data Validation 


and momentum flux ratio for two phase jet is given by 




(pl ul A fuel Pg U airbl ^ airbl 


spray 


p g u 


2 

cross 


The 120 micron droplet predicted the most penetration and the 80 micron droplet 
predicted the least penetration. 
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To validate the results obtained using the Excel spreadsheet (because no analytical 
solution is available), a FORTRAN program was also developed for this purpose and is 
given in the appendix along with three validation cases (plots). FORTRAN and Excel 
calculations are compared in Figure 5. With an Excel spreadsheet, we do not have to 
compile, build, and link as in a regular FORTRAN code. In addition, the graphics are 
immediately displayed after the computations are completed, so that the results are seen 
quickly and changes in the input can be made. The interactive spreadsheet is available on 
this CD as a separate document. Additional copies of the spreadsheet can be requested by 
e-mailing: Cecil.J.Marek@grc.nasa.gov . The report portion can be accessed on the web 
at: http://gltrs.grc.nasa.gov/cgi-bin/GLTRS/browse.pl72002/TM-2Q02-21 1710.html 
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User’s Manual 


This program is written in Microsoft Visual Basic Excel. There are three sheets in the 
program, namely the Instruction Sheet, the Process Sheet, and the Code Sheet. 

Instructions Sheet 


The instruction sheet contains a brief description of the problem, the solution method, 
and the user-input variables. Several schematics, which describe the forces applied on the 
droplet, can also be found on this sheet. 

Process Sheet 


The process sheet contains the user-inputs and the solution plots. There are five graphs on 
this sheet (see Figure 6), namely droplet trajectory, droplet velocity profile, drag 
coefficient, Cd, as a function of time, droplet velocity profiles as a function of time, and 
droplet trajectory profiles as a function of time. The cells highlighted in green are the user 
inputs. The cells highlighted in cyan contain computed values associated with the green 
cells; therefore, they are locked to prevent the user from modifying the values. When 
values are entered into the formula cells, the formulas are erased and linkages to other 
cells are interrupted. That is why for this version we chose to lock the closed cells 
without a password. 


Inputs 

Spherical radius of the droplet (^ 2 ) 
Viscosity of the crossflow fluid 
Gravitational Acceleration 
Droplet Density^ 

Fluid Density 
Time Step 

Fluid Velocity (X direction)'' 

Fluid Velocity (Z direction) 


Input Sheet 


0.00006 m 



Initial Droplet Posit ionjpi<^tfectiafi) X Q 

Initial Droplet VdHEity (X direanon) Ud Q 

Initial Droplet Position (Z dilution) Zo 

Initial Droplet Velocity (Z direction) Wd Q 


Projected area of th/droplet 
Volume of the drowet 
Total Cycles (IN^GER NUMBER) 
Data taken every (INTEGER NUMBER) 
Data plot for 


1.13ME-08 m 2 
9.0P4E-13 m 3 i 
Chiles 
300 Cycles 
6. OOE-O “^seconds; 


11314.29 /urn 2 
905142.9 ^tm 3 


Coordinate Systems 


X-RIGHT POSITIVE 
Z--UP POSITIVE 


COMMAND Update | Halt | Clear Data | 

COMPLETION 


i ► h\ CoverPage / Instruction \ Process /Code / 


Figure 5 
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Droplet Distance Profile 


Input Sheet 


Inputs 


Spherical radius of the droplet (r 32 ) 

r 32 = 

60 

m 

Viscosity of the crossflow fluid 

g= 

2.00E-05 

Kg/m-s 

Gravitational Acceleration 

g = 

9.8 

m/s 2 

Droplet Density 

d = 

822 

Kg/m 3 

Fluid Density 

g = 

1.22 

Kg/m 3 

Time Step 

t = 

1.00E-06 

sec 

Fluid Velocity (X direction) 

Ug = 

0 

m/s 

Fluid Velocity (Z direction) 

Wg = 

-38 

m/s 


or 0.00006 m 


3 

crS’ 

£ 

On 


Initial Droplet Position (X direction) X D = 

Initial Droplet Velocity (X direction) Ud 0 = 

Initial Droplet Position (Z direction) Z D = 

Initial Droplet Velocity (Z direction) Wd D = 

0 

-2.4 

0 

0 

m/s 

m 

m/s 




Projected area of the droplet A d = 

Volume of the droplet V d = 

Total Cycles (INTEGER NUMBER) 

Data taken every (INTEGER NUMBER) 

Data plot for 

1.13143E-08 

9.05143E-13 

m 2 or 11314.29 m 2 

m 3 or 905142.9 m 2 

600000 

300 

Cycles 

Cycles 

seconds; 2001 Data Written 

6.00E-01 


Coordinate Systems X-RIGHT POSITIVE 

Z-UP POSITIVE 


COMMAND 

COMPLETION 



OPTIONAL: For TecPIot Viewer ONLY 

Do you want to save the data into a file? 4 Yes 
3 No 


N 

3 


b 


1.2 

1 

0.8 

0.6 

0.4 

0.2 

0 


0 0.2 0.4 0.6 0.8 


Distance-X(m) 


Droplet Velocity Profile 



Velocity-X direction (m/s) 


Droplet Drag Coefficient C D Vs Time 


0.8 
0° 0.6 

0.4 

0.2 


0 


0 


0.2 0.4 0.6 0.8 

Time (sec) 








From Figure 5, the last two user inputs, Total Cycles (C25) and Data taken every ### 
cycles (C26) are included. The number assigned in cyan cell “E27”must be kept below 
65,536; the cell will turn into red if this condition is not satisfied. Cell E27 basically 
shows amount of data will be written into the code sheet, and Excel can only hold 65,536 
rows of data. Keeping the value below the limit can be done by changing the value in the 
cell “C26”. 

After all the inputs have been specified, clicking the Update button will instruct the 
program to update the data and the five plots. 

In summary, the user needs to do the following steps to run the program: 

1 . Go to the Process Sheet h by clicking on the Process tab 

2. Enter input values in the green cells 

3. Adjust the value in the cell C26 so that the computed value in cell E27 is less 
than 65,536 

4. Click the Update button 

5. Observe the droplet profiles on the five solution plots 

6. Repeat step 1 through step 5 for different input values 

7. Click the Clear Data button to clear the data (Optional). 

Additional features include the option to store the computed data into a TecPlot format 
file. This feature provides the user a flexibility to plot the data using other software such 
as TecPlot. The details about how to run this option are explained in the discussion 
section below. 

Code Sheet 


The last sheet is the code sheet, which contains the solution data produced by the 
program. Six columns namely A, B, C, D, E, and F hold the data of time, trajectory x- 
component, trajectory z-component, velocity x-component, velocity z-component, and 
drag coefficient respectively. 

Discussions 


The five plots are based on the X-Z coordinate system where X positive is to the right 
and Z positive is up. If the inputs of the droplet’s velocity injection is set only in the X 
positive component, and the crossflow’s velocities are set to zero, you will see a 
decreasing slope profile on the trajectory plot. Because relative velocity in the Z direction 
between the droplet and the crossflow is zero, gravity force overcomes zero buoyancy 
force, and results in the droplet to move down, in the negative Z direction. 

The total number of cycles will determine the time for the program to process the data. 
Data taken every cycle also determines the speed of the program to process the data. The 
more data is collected, the slower the program to complete the cycle. 
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You are encouraged to change the chart type of the five plots which you find best and 
easy to analyze. You might want to switch to dotted points-connected line chart to 
analyze the data points. 

The cells highlighted in cyan and in white are locked for the coding security purpose. It is 
highly recommended that initially the user not unlock and modify these cells. Any 
modifications made may cause the program to crash. 

The last feature added to this program is the flexibility for the user to store the computed 
data into TecPlot format file. By enabling this option, the user needs to specify the path 
and the filename to store the output. This option must be disabled if the path and the file 
name are not specified. 
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Appendix A 
Visual Basic Code 


Under the ‘Tools’, ‘Macro’, ‘Visual Basic’ button, the heart of the numerical code is 
presented. This code is reproduced here in case something happens to the program. Its 
logic is similar to the FORTRAN code following, but some things are different. 


Private Sub CommandButton3_Click() 

Dim n As Long, nn As Double, nnnn As Long, userchoice As Long 
Dim xx(0 To 1) As Double, xp(0 To 1) As Double 
Dim zz(0 To 1) As Double, zp(0 To 1) As Double 
Dim Cdm As Double, Rems As Double 

Dim A_x As Double, A_z As Double, h As Double, kl As Double, r As Double, ht As 
Double 

Dim kl As Double, k2 As Double, k3 As Double, k4 As Double, 11 As Double, 12 As 
Double, 13 As Double, 14 As Double 

Dim kzl As Double, kz2 As Double, kz3 As Double, kz4 As Double, lzl As Double, 
lz2 As Double, lz3 As Double, lz4 As Double 
Dim Urm As Double 

Dim ug As Double, wg As Double, mu As Double, rg As Double 
Dim a As Double, b As Double 
Dim location 
Dim sd, sf 

Dim al, bl, cl, dl, el, fl 
Dim unitconv 

Call Macro2 
Halt = False 

Msg = "Do you want to continue ?" 

Style = vbYesNo 

Title = "Jet Flow in CrossFlow" 

'Getting the input 

ug = Sheets("Process").Cells(l 1, 3) 
wg = Sheets("Process").Cells(12, 3) 
mu = Sheets("Process").Cells(6, 3) 
rg = Sheets("Process").Cells(9, 3) 
rd = Sheets("Process").Cells(8, 3) 
g = Sheets("Process").Cells(7, 3) 
xx(0) = Sheets("Process").Cells(17, 3) 
xp(0) = Sheets("Process").Cells(18, 3) 
zz(0) = Sheets("Process").Cells(19, 3) 
zp(0) = Sheets("Process").Cells(20, 3) 
h = Sheets("Process").Cells(10, 3) 
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kl = Sheets("Process").Cells(25, 3) 
userchoice = Sheets("Process").Cells(26, 3) 
r = Sheets("Process").Cells(5, 3) 
r = r/ 1000000# 


Worksheets("Process").CommandButtonl. Width = 0 
Worksheets("Process").CommandButtonl. Visible = True 

n = 0 
nn = -1 
nnn = 1 
ht = 0 


20 

ht = ht + h 
a = xp(n) 
b = zp(n) 

Call Runge(a, b, r, rd, rg, g, mu, ug, wg, A_x, A_z) 
If (SuperExit) Then GoTo 50 
kl = h * xp(n) 

11 = h* A_x 
kzl = h * zp(n) 
lzl = h * A_z 

a = xp(n) + 11/2# 
b = zp(n) + lzl / 2# 

Call Runge(a, b, r, rd, rg, g, mu, ug, wg, A_x, A_z) 
If (SuperExit) Then GoTo 50 
k2 = h * xp(n) 

12 = h * A_x 
kz2 = h * zp(n) 
lz2 = h * A_z 

a = xp(n) + 12/2# 
b = zp(n) + lz2 / 2# 

Call Runge(a, b, r, rd, rg, g, mu, ug, wg, A_x, A_z) 
If (SuperExit) Then GoTo 50 
k3 = h * xp(n) 

13 = h * A_x 
kz3 = h * zp(n) 
lz3 = h * A_z 

a = xp(n) + 13 
b = zp(n) + lz3 

Call Runge(a, b, r, rd, rg, g, mu, ug, wg, A_x, A_z) 
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If (SuperExit) Then GoTo 50 
k4 = h * xp(n) 

14 = h * A_x 
kz4 = h * zp(n) 
lz4 = h * A_z 

xp(n + 1) = xp(n) + (1# / 6#) * (11 + 2# * 12 + 2# * 13 + 14) 
zp(n + 1) = zp(n) + (1# / 6#) * (lzl + 2# * lz2 + 2 #* lz3 + lz4) 

xx(n + 1) = xx(n) + (1# / 6#) * (kl + 2# * k2 + 2# * k3 + k4) 
zz(n + 1) = zz(n) + (1# / 6#) * (kzl + 2# * kz2 + 2# * kz3 + kz4) 

nn = nn + 1# 

If nn = 0 Then 

Urm = (((xp(n) - ug) A 2) + ((zp(n) - wg) A 2)) A 0.5 
Rems = 2 * rg * Urm * r / mu 
If Rems <= 1000 Then 

Cdm = (((Rems A (2 / 3)) / 6) + 1) * 24 / Rems 
Else 

Cdm = 0.424 
End If 

Sheets("Code").Cells(2, 1) = nn * h 
Sheets("Code").Cells(2, 2) = xx(n) 

Sheets("Code").Cells(2, 3) = zz(n) 

Sheets("Code").Cells(2, 4) = xp(n) 

Sheets("Code").Cells(2, 5) = zp(n) 

Sheets("Code").Cells(2, 6) = Cdm 

Worksheets("Process").CommandButtonl.Width = ((nn/kl) * 165.75) 
Worksheets("Process").CommandButtonl. Caption = ((nn / kl) * 100) & "%" 
Worksheets("Process").CommandButtonl. Height = 20.25 

Else 

If 0 = (nn Mod userchoice) Then 
nnn = nnn + 1 

Urm = (((xp(n + 1) - ug) A 2) + ((zp(n + 1) - wg) A 2)) A 0.5 
Rems = 2 * rg * Urm * r / mu 
If Rems <= 1000 Then 

Cdm = (((Rems A (2 / 3)) / 6) + 1) * 24 / Rems 
Else 

Cdm = 0.424 
End If 

Sheets("Code").Cells(nnn + 1, 1) = nn * h 
Sheets("Code").Cells(nnn +1,2) = xx(n + 1) 

Sheets("Code").Cells(nnn +1,3) = zz(n +1) 

Sheets("Code").Cells(nnn + 1,4) = xp(n + 1) 
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Sheets("Code").Cells(nnn +1,5) = zp(n +1) 

Sheets("Code").Cells(nnn +1,6) = Cdm 

Worksheets("Process").CommandButtonl.Width = ((nn/kl) * 165.75) 
Worksheets("Process").CommandButtonl. Caption = ((nn / kl) * 100) & "%" 
Worksheets("Process").CommandButtonl. Height = 20.25 

DoEvents 
If (Halt) Then 
DoEvents 
Halt = False 

Response = MsgBox(Msg, Style, Title) 

If Response = vbNo Then 
GoTo 10 
Else 
End If 
End If 

Else 
End If 
End If 


xp(n) = xp(n + 1) 
zp(n) = zp(n +1) 
xx(n) = xx(n + 1) 
zz(n) = zz(n +1) 

If nn >= kl Then 
GoTo 10 
Else 

GoTo 20 
End If 


10 

Macro 1 

sd = Sheets("Process").Cells(27, 5) 

If Optional_Save = "YES" Then 
location = Sheets("Process").Cells(44, 2) 

Open location For Output As #1 

Print #1, "TITLE = "; Spc(2); ; "Spray Jet In CrossFlow"; 

If Unit = "millimeter" Then 
unitconv = 1000# 

Print #1, "Variables ="; Spc(2); ; "Time(sec)"; ; Spc(2); ; "X(mm) M ; "" 

Spc(2); ; "Z(mm)"; Spc(2); ; "U(mm/s) M ; ; Spc(2); ; "W(mm/s)" 

; Spc(2); ; "Cd"; 

Else 
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unitconv =1# 

Print #1, "Variables ="; Spc(2); "Time(sec)"; ; Spc(2); ; "X(meter)"; 

Spc(2); ; "Z(meter)"; ; Spc(2); ; "U(m/s) M ; ; Spc(2); ; ”W(m/s)"; 

Spc(2); ; "Cd"; 


9 

ft ft ft ft . 


End If 

Print #1, "ZONE I ="; sd & Spc(2); "F = POINT" 


For sf = 2 To sd + 1 
al = Sheets("Code").Cells(sf, 1) 
bl = Sheets("Code").Cells(sf, 2) 
bl = bl * unitconv 
cl = Sheets("Code").Cells(sf, 3) 
cl = cl * unitconv 
dl = Sheets("Code").Cells(sf, 4) 
dl = dl * unitconv 
el = Sheets("Code").Cells(sf, 5) 
el = el * unitconv 
fl = Sheets("Code").Cells(sf, 6) 

Print #1, al; Spc(2); bl; Spc(2); cl; Spc(2); dl; Spc(2); el; Spc(2); fl 
Next sf 
Close #1 
Else 


End If 
50 


SuperExit = False 
End Sub 


Sub Runge(xps, zps, rs, rds, rgs, g, mus, ugs, wgs, A_xs, A_zs) 

Dim Area As Double, Volume As Double, Ur As Double, Re As Double, Cd As Double 

Area = (22# / 7#) * rs A 2 

Volume = (22# / 7#) * (4# / 3#) * rs A 3 

Ur = (((xps - ugs) A 2#) + ((zps - wgs) A 2#)) A 0.5 

Msg = "The Time Step Size is too big. Program Terminated Ur = " & Ur 

Msg2 = "Reduce the time step size!" 

Style = vbOKOnly 

Title = "Jet Flow in CrossFlow" 

If Ur > (10 A 80) Then 

Response = MsgBox(Msg, Style, Title) 

Response = MsgBox(Msg2, Style, Title) 

SuperExit = True 

Sheets("Process").Range("C 1 0"). Select 
Else: End If 

Re = 2# * rgs * Ur * rs / mus 
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If Re <= 1000 Then 
Cd = Re A (2# / 3#) 

Cd = Cd / 6# 

Cd = (Cd + 1#) * 24# / Re 
Else 


Cd = 0.424 
End If 

A_xs = -0.5 * rgs * Ur * Area * Cd * (xps - ugs) / (rds * Volume) 

A_zs = (-0.5 * rgs * Ur * Area * Cd * (zps - wgs) + (rgs - rds) * Volume * g) / (rds * 
Volume) 

End Sub 
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Appendix B 

Equivalent Fortran Program 


*deck main 


c This program calculates droplet's trajectory and its velocity * 
c after it is injected into a cross-stream using the 4th-order * 
c Runge-Kutta method. It calls derivs which defines the ordinary * 
c differential equations (ODEs) . * 

c All the constant values can be found in subroutine derivs * 


c 

parameter (nd=10000000 ) 

dimension x (nd) , dx (nd) , u (nd) , du (nd) , t (nd) , xt (nd) , ut (nd) , 

& dxt (nd) , dut (nd) , dxm (nd) , dum (nd) 

dimension z (nd) , dz (nd) , w (nd) , dw (nd) , zt (nd) , wt (nd) , 

Sc dzt (nd) , dwt (nd) , dzm (nd) , dwm (nd) 

c + + + 

c +++ input parameters 
c + + + 

print*, 'Enter time step : ' 
read* , dt 
print* , ' ' 

print*, 'Enter total time: ' 
read* , tend 

print*, 'Enter initial droplet position x(0) : ' 

read* , x (1) 
print* , ' ' 

print*, 'Enter initial droplet velocity u(0) : ' 

read* , u ( 1 ) 
print* , ' ' 

print* ,' Enter initial droplet position z(0) :' 

read* , z (1) 
print* , ' ' 

print*, 'Enter initial droplet velocity w(0) : ' 

read* , w ( 1 ) 
print* , ' ' 

print*, 'Enter radius [micro-meter] : ' 
read* , r 
print* , ' ' 

r=r* 1 . e-6 
c 

pi=acos ( - 1 . ) 

ad=pi*r**2 ! Droplet area 

vd=4 . /3 . *pi*r**3 ! Droplet volumn 

c 

print*, 'Area = ' , ad 
print*, 'Volumn = ' ,vd 
c 

t (1) =0 . 0 
c 

n=tend/ dt 
c 
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write ( * , * ) 

write (*,*) ' time step dt 
write (* , * ) 


, dt 


if (n . gt . nd) then 

write (*,*) 1 n is greater than nd -- run aborting' 
write (*,*)' n = ' ,n, ' nd = ' , nd 
stop 
endif 
c 

t (n) =tend 
c + + + 

c +++ Solving the differential equations using 4th-order 
c +++ Runge Kutta methods, 
c + + + 

dtt=0 . 5*dt 
do 10 i=l , n- 1 
c + + + 
c + + + 

call derivs ( i , r , ad, vd, u, dx, du, w, dz , dw) 
c + + + 
c + + + 

xt ( i ) =x ( i ) +dtt*dx ( i ) 
ut ( i ) =u ( i ) +dtt*du ( i ) 
zt ( i ) =z ( i ) +dtt*dz ( i ) 
wt ( i ) =w ( i ) +dtt*dw ( i ) 

call derivs (i , r , ad, vd, ut , dxt , dut , wt , dzt , dwt ) 
c + + + 
c + + + 

xt ( i ) =x ( i ) +dtt*dxt ( i ) 
ut ( i ) =u ( i ) +dtt*dut ( i ) 
zt ( i ) =z ( i ) +dtt*dzt ( i ) 
wt ( i ) =w ( i ) +dtt*dwt ( i ) 

call derivs ( i , r , ad, vd, ut , dxm, dum, wt , dzm, dwm) 
c + + + 
c + + + 

xt ( i ) =x ( i ) +dt*dxm ( i ) 
ut ( i ) =u ( i ) +dt*dum ( i ) 
dxm ( i ) =dxt ( i ) +dxm ( i ) 
dum ( i ) =dut ( i ) +dum ( i ) 
zt ( i ) =z ( i ) +dt*dzm ( i ) 
wt ( i ) =w ( i ) +dt*dwm ( i ) 
dzm ( i ) =dzt ( i ) +dzm ( i ) 
dwm ( i ) =dwt ( i ) +dwm ( i ) 

call derivs (i , r , ad, vd, ut , dxt , dut , wt , dzt , dwt) 
c + + + 

c +++ advance the solutions to n+1 time step 
c + + + 

t ( i+1 ) =t ( i ) +dt 

x ( i+1 ) =x ( i ) +dt/ 6 . * (dx ( i ) +dxt ( i ) +2 . *dxm ( i ) ) 
u ( i + 1 ) =u ( i ) +dt/ 6 . * (du ( i ) +dut ( i ) +2 . *dum ( i ) ) 
z ( i+1 ) =z ( i ) +dt/ 6 . * (dz ( i ) +dzt ( i ) +2 . *dzm ( i ) ) 
w ( i+1 ) =w ( i ) +dt / 6 . * (dw ( i ) +dwt ( i ) +2 . *dwm ( i ) ) 
c 

10 continue 
c 

c + + + 
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c +++ write data sets to file result.dat 
c + + + 

open (unit=3 , f ile= ' result . dat ' , f orm= ' formatted ' ) 
do 20 i=l,n 

write (3,100) t(i),x(i),z(i),u(i),w(i) 

20 continue 
c 

write (*,*) ' The data is written to result.dat as' 
write(*,*) ' t, x, z, u, w' 
c 

close (3 ) 
c 

100 format ( lx, f 5 . 3 , 4 ( lx, f 15 . 10 ) ) 
end 
c 

c=================================================== 

*deck derivs 

subroutine derivs ( i , r , ad, vd, u, dx, du, w, dz , dw) 


c This subroutine defines the differential equations 
c being solved by the program. * 


dimension u ( 1 ) , dx ( 1 ) , du ( 1 ) , w ( 1 ) , dz ( 1 ) , dw ( 1 ) 

c + + + 

c +++ Constant values 
c + + + 

rhod=300. ! Droplet density 

rhog=1.19 ! Fluid density 

ug=0 . ! Fluid velocity in x-direction 

g=9 . 8 ! Gravity 

wg=-5.0 ! Fluid velocity in z-direction 

eni=0.0004 ! viscosity of the gaseous fluid 

c 

ur=abs (sqrt ( (u(i) -ug) **2 + (w(i) -wg) **2) ) 
c 

red=2 . *r*ur*rhog/eni 
if (red . le . 1000 . ) then 
cd=24 . / red* ( 1 . +1 . /6 . *red* * (2./3.) ) 
else 

cd=0 .424 
endif 
c 

dx ( i ) =u ( i ) 

du ( i ) =1 . / (rhod*vd) * ( - 0 . 5*rhog* (u ( i ) -ug) *ur*ad*cd) 
c 

dz ( i ) =w ( i ) 

dw ( i ) =1 . / (rhod*vd) * ( - 0 . 5*rhog* (w ( i ) -wg) *ur*ad*cd+ 
& (rhog-rhod) *vd*g) 
c 

return 

end 
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